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ABSTRACT 

We report the discovery of HAT-P-16b, a transiting extrasolar planet orbiting the V = 10.8 mag F8 
dwarf GSC 2792-01700, with a period P = 2.775960 ± 0.000003 d, transit epoch T c = 2455027.59293 ± 
0.00031 (BJD***), and transit duration 0.1276±0.0013d. The host star has a mass of 1.22 ±0.04 M , 
radius of 1.24±0.05i? Q , effective temperature 6158±80K, and metallicity [Fe/H] = +0.17±0.08. The 
planetary companion has a mass of 4.193 ± 0.094 Mj, and radius of 1.289 ± 0.066 Rj yielding a mean 
density of 2.42±0.35gcm -3 . Comparing these observed characteristics with recent theoretical models, 
we find that HAT-P-16b is consistent with a 1 Gyr H/He-dominated gas giant planet. HAT-P-16b 
resides in a sparsely populated region of the mass-radius diagram and has a non-zero eccentricity of 
e = 0.036 with a significance of 10<r. 

Subject headings: planetary systems — stars: individual (HAT-P-16, GSC 2792-01700) techniques: 
spectroscopic, photometric 



1. INTRODUCTION 

Planets that transit their host stars play a special role 
in our understanding of the characteristics of exoplanets: 
their transit allows us to accurately determine the radius 
and the orbital inclination of the planet from the photo- 
metric light curve so that an actual mass can be derived 
from a spectroscopic orbit of the host star. The mass and 
radius enable us to infer a bulk composition of the planet, 
and although there are degeneracies associated with the 
bulk composition, it allows us to put constraints on mod- 
els of planetary structure and formation theories. The 
incredible diversity of the over 60 discovered transiting 
planets, ranging from dense planets with a higher mean 
density than that of copper to strongly irradiated puffed- 
up planets with a mean density comparable to that of 
corkwood, have baffled the exoplanet community, and 
no unified theory has been established to explain all the 
systems consistently. Transiting extrasolar planet (TEP) 
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discoveries are primarily the result of ded icated ground- 
based searches, such as SuperW ASP (Pollacc o et al 
2006), HATNet dBakos et al.ll200l TrES (jAlonso et al 
■20041 ) and XO (McCullough ct al. 200 51), and sp ace- 
borne s earches, such as CoR oT (jBaglin et al.ll2006| ) and 
Kepler (TBorucki et al.H201Ct ). 

Since its commissioning in 2003, the Hun garian-made 
Auto mated Telescope Network (HATNet; iBakos et al.l 
|2004|) survey has been one of the major contributors 
to the discoveries of TEPs. HATNet has discovered 
over a dozen TEPs since 2006 by surveying bright stars 
(8 ;$ I ?$ 12.5 mag) in the Northern hemisphere and 
has now covered approximately 11% of the Northern sky. 
HATNet consists of six wide field automated telescopes; 
four of these are located at the Fred Lawrence Whipple 
Observatory (FLWO) in Arizona, and two on the roof 
of the Submillimeter Array hangar (SMA) of the Smith- 
sonian Astrophysical Observatory (SAO) in Hawaii. In 
this paper we report a new TEP discovery of HATNet, 
called HAT-P-16b. 

The structure of the paper is as follows. In Section [2] 
we summarize the observations, including the photomet- 
ric detection, and follow-up observations. In Section [3] 
we describe the analysis of the data, such as the stellar 
parameter determination (Section 13. 1[) . ruling out blend 
scenarios (Section 13.21) . and global modeling of the data 
(Section [33)) . We discuss our findings in Section [4] 

2. OBSERVATIONS 

2.1. Photometric detection 

The transits of HAT-P-16b were detected with the 
HAT-6 and HAT-7 telescopes in Arizona and the HAT-8 
and HAT-9 telescopes in Hawaii. The star GSC 2792- 
01700 lies in the intersection of 3 separate HATNet fields 
internally labeled as 123, 163 and 164. Field 123 was ob- 
served with an / filter and a 2Kx2K CCD, while fields 
163 and 164 were observed through R filters with 4Kx4K 
CCDs. All three fields were observed with a 5 minute ex- 
posure time and at a 5.5 minute cadence in the period be- 
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Fig. 1.— Unbinned light curve of HAT-P-16 including all 12,552 
instrumental / band and R band 5.5 minute cadence measure- 
ments obtained with the HAT-6, HAT-7, HAT-8 and HAT-9 tele- 
scopes (see the text for details), and folded with the period of 
P = 2.7759602 days (which is the result of the fit described in 
Section [3}. The / band and R band data are merged in this fig- 
ure. Th e sol id line shows the transit model fit to the light curve 
(Section 13.31 . In the lower panel, the large solid circles represent 
the binned average of the photometric measurements with a bin 
size of 0.005 days. 

tween July 2007 and September 2008 during which over 
12,000 exposures were gathered for the three fields. Each 
image contains between 27,000 and 76,000 stars down to 
a magnitude of / ~ 13 for field 123 and R ~ 15 for 
fields 163 and 164, yielding a photometric precision for 
the brightest stars in the field of ~ 2 mmag for field 123 
and ~ 5 mmag for fields 163 and 164. 

Frame calibration, astrometry and aperture photom- 
etry were done in an identical w ay to recent HAT - 
Net TEP disco veries, as described in feakos et all (2009) 
and Pal (2009). The resulting light curves were decor- 
related (cleaned of trends) using the External Pa- 
rameter De correlation techn ique in "constant" mode 
(EPD; see iBakos et all 12009ft and the Trend Filter- 
ing Algorithm (TFA; see IKovacs et al.l 120051 ). The 
light curves were searched for periodic box-like sig- 
nals using the Bo x Least-Squares method (BLS; see 
IKovacs et al.l [2002) . A significant signal was detected 
in the light curve of GSC 2792-01700 (also known as 
2MASS 00381756+4227470; a = 00 h 3 8 m 17.56s 5 = 
+42°27'47.1"; J2000; V = 10.812: iDroege et al.ll200fih . 
with a depth of ~ 10.2 mmag, and a period of P = 
2.7760 days. The dip had a relative duration (first to 
last contact) of q rj 0.0460 ± 0.0005, corresponding to a 
total duration of Pq w 3.062 ± 0.031 hr (see Figured]). 

2.2. Reconnaissance Spectroscopy 

Transiting planet candidates found by ground-based, 
wide-field photometric surveys must undergo a rigorous 
vetting process to eliminate the many astrophysical sys- 
tems mimicking transiting planets (called false positives), 
the rate of which has proved to be much higher than the 
occurrence of true planets (10 to 20 times higher). Low 
signal to noise ratio (S/N) high-resolution reconnaissance 
spectra are used to extract stellar parameters such as ef- 
fective temperature, gravity, metallicity and rotational 
and radial velocities to rule out these false positives. Ex- 
amples of false positives which are discarded are eclipsing 
binaries and triple systems. The latter can be either hi- 
erarchical or chance alignment systems where the light 
of the eclipsing pair of stars is diluted by the light of 



a third brighter star. Rapidly rotating and/or hot host 
stars whose spectrum is unsuitable for high precision ve- 
locity work are also discarded. 

We acquired 7 reconnais sance spectra with the CfA 
Digital Speedometer (DS; Latham 1992) mounted on 
the FLWO 1.5 m Tillinghast Reflector between De- 
cember 2008 and January 2009. The extracted 
modest precision radial velocities gave a mean abso- 
lute RV=— 16.83kms _1 with an rms of 0.51 kms -1 , 
which is consistent with no detectable RV varia- 
tion. The stellar parameters derived from the spectra 
([Torres. Neuhauser. fe; Guen thcr 2002), including the ef- 
fective temperature Teg* = 6000 ± 100 K, surface gravity 
log g+ = 4.0 ± 0.25 (log cgs) and projected rotational 
velocity vsmi = 3.8 ± 1.0 km s -1 , correspond to a F8 
dwarf. 

2.3. High resolution, high S/N spectroscopy 

The high significance of the transit detection by HAT- 
Net, together with the stellar spectral type and small 
RV variations encouraged us to gather high-resolution, 
high S/N spectra to determine the orbit of the sys- 
tem. We have taken 21 spectra between August and 
October 2009 using the Fibre-fed Echelle Spectrograph 
(FIES) at the 2.5 m Nordic Optical Telescope (NOT) 
at La Palma, Spain (Djupvik fc Anders en 2009). We 
used the medium and the high-resolution fibers (l'/3 pro- 
jected diameter) with resolving powers of A/AA « 46,000 
and 67,000, respectively, giving a wavelength coverage of 
~ 3600 - 7400 A. Recently, FIES has also been used to 
obtain a spectroscopic orbi t for one of the firs t Kepler 
planets, namely Kepler- 7b (Lat ham et al.ll2010D. 

We also used the HIRES instrument (|Vogt et al.lll994ft 
on the Keck I telescope located on Mauna Kea, Hawaii. 
We obtained 6 exposures with an iodine cell, plus 1 
iodine-free template with Keck-I/HIRES. The observa- 
tions were made on 6 nights between 2009 July 3 and 
2009 October 29. The width of the spectrometer slit 
used on HIRES was 0'.'86, resulting in a resolving power 
of A/AA w 55,000, with a wavelength coverage of ~ 
3800 — 8000 A. The iodine gas absorption cell was used 
to superimpose a dense forest of I2 lines on the stellar 
spec trum and establish an accurate wavelength fiducial 
(see iMarcv fc Butlerlfl992l ). Relative RVs in the Solar 
System bary centric frame were derived as described by 
iButler et al.l ijl996), incorporating full modeling of the 
spatial and temporal variations of the instrumental pro- 
file. 

The final RV data and their errors (for both instru- 
ments) are listed in Table Q] The folded data are plotted 
in Figured The systemic gamma velocities (that are the 
result of the global modeling, as laid out in Section [3]) 
have been subtracted to ensure a common zero-point. 
The best orbital fit (see Section |3|) is superimposed in 
the figure. 

Since this is the first time the NOT has been used 
to determine the orbit of a HATNet planet, we describe 
briefly the extraction of the radial velocities. We have 
built a custom pipeline to rectify and cross correlate 
spectra from Echelle spectrographs with the goal of pro- 
viding very precise radial velocities. The first step in 
the extraction process is to remove the bias level and 
crop the raw images. To effectively remove cosmic rays, 
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Fig. 2. — Top: RV measurements from the Nordic Optical Tele- 
scope and Keck for HAT-P-16, along with an orbital fit, shown as a 
function of orbital phase, using our best-fit period (see Section Op- 
Velocities from NOT using the medium and high resolution fiber 
are shown as squares and triangles, respectively, and the Keck ve- 
locities are shown as circles. Zero phase is defined by the transit 
midpoint. The center-of-mass velocity has been subtracted. Note 
that the error bars include the stellar jitter, added in quadrature to 
the formal errors given by the spectrum reduction pipeline. Second 
panel: phased residuals after subtracting the orbital fit (also see 
Section[3]l. The rms variation of the residuals is about 10.0ms -1 . 
Third panel: bisec tor s pan variations (BS) including the template 
spectrum (Section l3.2l l. Bottom: relative S activity index values 
which have only been calculated for the Keck spectra. Note the 
different vertical scales of the panels. 



the observation is divided into three separate exposures, 
which enables us to combine the raw images using me- 
dian filtering, removing virtually all cosmic rays. We 

use a flat field to trace the Echelle orders and we rectify 
the spectra by using the "optima l extraction" algorithm 
(jHewett et all 1 198a IHorndfl9860 . The blaze function is 
determined by fitting cubic splines to a combined high 
S/N flat field exposure and is saved separately in order 
to preserve the original flux in the stellar exposure. By 
dividing the normalized blaze function into the rectified 
flat field spectrum, we can determine the pixel to pixel 
variations and correct for these. The scattered light in 
the two-dimensional raw image is determined and re- 
moved by masking out the signal in the Echelle orders 
and fitting the inter-order scattered light flux with a 
two-dimensional polynomial. 

Thorium argon (ThAr) calibration images are obtained 
through the science fiber before and after the stellar ob- 
servation. The two calibration images are combined to 
form the basis for the fiducial wavelength calibration. 
We have determined that the best wavelength solution 



is achieved by choosing an exposure time that saturates 
the strong argon lines but enhances the forest of weaker 
thorium lines. 

FIES is not a vacuum spectrograph and is only tem- 
perature controlled to 0.1 °C. Consequently, the radial 
velocity errors are dominated by our limited ability to 
correct for shifts due to pressure, humidity and temper- 
ature variations. In order to successfully remove these 
large variations (> 1.5kms -1 ), it is critical that the 
ThAr light travels through the same light path as the 
stellar light and thus acts as an effective proxy to remove 
these variations. We have therefore chosen to interleave 
the stellar observation between two ThAr exposures in- 
stead of using the simultaneous ThAr technique, which 
may not exactly describe the induced shifts in the science 
fiber. The centers of the ThAr lines are found by fitting a 
Gaussian function to the line profiles and a two- dimen- 
sional fifth order Legendre polynomial is used to describe 
the wavelength solution. 

Once the spectra have been extracted, a multi-order 
cross correlation is performed order by order. First, the 
spectra are interpolated to a common oversampled log 
wavelength scale with the same monotonic wavelength 
increment in all orders. A high and low pass filter is ap- 
plied in the Fourier domain and the ends of the spectra 
are apodized with a cosine bell function. The orders are 
cross correlated using a Fast Fourier Transform (FFT) 
and the cross correlation function (CCF) for each order 
is co-added. This automatically weights each order by 
its flux, giving more weight to the orders with more pho- 
tons. This CCF peak is fitted with a Gaussian function 
to determine its center. The internal precision is esti- 
mated by finding the radial velocity for each order and 
the RV error is thus a = RMS(v)/vN, where v is the 
RV of the individual orders and N is the number of or- 
ders. For the final radial velocities, a template spectrum 
is constructed by shifting and co-adding all the observed 
spectra, and the individual spectra are cross correlated 
against this co-added template spectrum to minimize the 
contribution of noise in the template. 

2.4. Photometric follow-up observations 

To confirm the transit signal and obtain high-precision 
light curves for modeling the system, we conducted 
photometric follow-up observations with the KeplerCam 
CCD on the FLWO 1.2 m telescope. We observed 4 tran- 
sit events of HAT-P-16 on the nights of 2009 Septem- 
ber 10, September 21, October 19 and October 30 (Fig- 
ure |3J). On the four nights, 303, 181, 293 and 471 frames 
were acquired with a cadence of 41, 52, 41 and 32 seconds 
(30, 40, 30 and 20 seconds of exposure time) in the Sloan 
i band, respectively. 

The reduction of the images, including frame calibra- 
tion, astromet ry, and p hotometry were performed as de- 
scribed in iBakos et al.l (2009). We also performed EPD 
and TEA to remove trends simultaneously with the light 
curve modeling (for more details, see Section [3]). The fi- 
nal light curves are shown in the upper plots of Figure [31 
superimposed with our best-fit transit light curve mod- 
els (see also Section [3]); the photometry is provided in 
Tabled 

3. ANALYSIS 
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TABLE 1 
Relative radial velocity and bisector span variation measurements 

of HAT-P-16. 



BJD RV a 


CTRV 


BS 


°F3S 


Inst. 


(2,455,000+) (ms" 1 ) 


(ms" 1 ) 


(ms" 1 ) 


(ms— 1 ) 




048.64476 . 


-137.7 


4.6 


-13.4 


10.9 


FIES h 


049.70638 .. 




-339.5 


10.5 


7.5 


17.5 


FIES h 


050.68239 .. 




-911.9 


6.1 


24.1 


13.6 


FIES h 


051.73142 .. 




60.3 


5.8 


-0.7 


10.4 


FIES h 


052.68956 .. 




-579.8 


4.4 


8.7 


8.9 


FIES h 


053.71751 .. 




-674.7 


11.3 


-6.1 


22.8 


FIES h 


054.71713 .. 




66.2 


5.7 


-0.9 


8.0 


FIES h 


107.58955 .. 




0.7 


5.2 


-16.8 


10.1 


FIES m 


108.62791 .. 




-965.2 


5.3 


-1.2 


9.2 


FIES m 


109.56094 .. 




-310.4 


7.0 


-7.3 


16.0 


FIES m 


110.54715 .. 




-128.9 


7.5 


8.3 


8.1 


FIES m 


111.56509 .. 




-989.4 


4.1 


-9.8 


5.7 


FIES m 


112.60304 .. 




-53.7 


5.5 


-8.7 


10.4 


FIES m 


113.50725 .. 




-294.7 


1.7 


-27.4 


8.6 


FIES m 


114.53303 .. 




-924.1 


4.5 


-8.1 


10.7 


FIES m 


116.52227 .. 




-558.0 


5.3 


-8.0 


10.8 


FIES m 


122.51752 .. 




-941.1 


7.5 


-5.2 


11.8 


FIES m 


123.46870 .. 




-262.4 


1.7 


-9.3 


10.2 


FIES m 


124.48749 .. 




-161.9 


7.9 


-7.2 


14.7 


FIES m 


125.44928 .. 




1008.8 


5.9 


4.0 


11.9 


FIES m 


126.47958 .. 




-29.8 


5.2 


-10.1 


9.4 


FIES m 


017.09285 .. 




-531.7 


2.3 


-3.1 


4.3 


HIRES 


019.09413 .. 




179.5 


2.6 


-18.1 


5.8 


HIRES 


107.08869 .. 








-10.3 


1.6 


HIRES 


107.09645 .. 




445.3 


2.1 


-2.9 


3.9 


HIRES 


108.97959 .. 




-497.6 


2.6 


33.5 


1.7 


HIRES 


112.12355 .. 




-109.3 


2.8 


6.7 


2.1 


HIRES 


135.02181 .. 




510.9 


2.5 


-5.8 


1.6 


HIRES 




2009 Sep 10 
%&t 



2009 Sep 21 



2009 Oct 19 



.". _ 2009 Oct 30 

■'iV.*. 

■£*^ 



. i r** **. . . • 






Note. — For the iodine- free template exposure we do not measure the RV 
but do measure the BS. Such template exposures can be distinguished by the 
missing RV value. The Inst, column refers to the instrument used, i.e. the 
FIES spectrograph at the NOT using the medium and high resolution fibers 
or the HIRES spectrograph at Keck I. orv and ctbs are formal statistical 
errors. 

The fitted zero-point that is on an arbitrary scale (denoted as -y re ] EiG. 3. — 




-0.15 



-0.1 



0.1 



0.15 



Section |3- 3D has not been subtracted from the velocities. 



TABLE 2 
Follow-up photometry of HAT-P-16 



BJD 




Mag a 


"Mag 


Mag(orig) b 


Filter 


(2,400,000 


+) 










55085.69759 




0.00368 


0.00082 


9.93719 


i 




55085.69808 




0.00111 


0.00082 


9.93426 






55085.69860 




-0.00132 


0.00082 


9.93260 


j 




55085.69911 




-0.00200 


0.00082 


9.93017 






55085.69962 




0.00022 


0.00082 


9.93280 


i 



Note. — This table is available in a machine-readable form in 
the online journal. A portion is shown here for guidance regarding 
its form and content. 

a The out-of-transit level has been subtracted. These magnitudes 
have been subjected to the EPD and TFA procedures, carried out 
simultaneously with the transit fit. 

b These are raw magnitude values without application of the EPD 
and TFA procedures. 



3.1. Properties of the parent star 

The stellar atmospheric parameters were derived from 
the template spectrum obtained with the Keck/HIRES 
instrument. We used th e Spectroscopy Made Eas y 
(SME) analysis package of iValenti &; Piskunovl fl 996). 
along with the atomic-line database of IValenti fc Fischerl 
(2005) . This yielded the following initial values and un- 



-0.05 0.05 

Time from transit center (days) 

Unbinned instrumental i band transit light curves, 
acquired with KeplerCam at the FLWO 1.2 m telescope. Superim- 
posed are the best-fit transit model light curves. At the bottom of 
the figure we show the residuals from the fit. Error bars represent 
the photon and background shot noise, plus the readout noise. 

certainties (which we have conservatively increased to 
include our estimates of the systematic errors): effec- 
tive temperature T c s± — 6175 ± 80 K, stellar surface 
gravity \ogg± = 4.44 ± 0.10 (cgs), metallicity [Fe/H] = 
0.15±0.06dex, and projected rotational velocity vsini = 
4.4±0.5kms" 1 . 

We could now use the effective temperature and the 
surface gravity found by SME to determine other stel- 
lar parameters, such as the mass, M*, and radius, i?*, 
using model isochrones. However, there may be strong 
correlations between effective temperature, gravity and 
metallicity in the spectroscopic determination of these 
parameters. Also, the effect of log g+ on the spectral line 
shapes is typically subtle, and as a result log g+ is gener- 
ally a rather poor luminosity indicator. Instead, we used 
the a/R+, the normalized semimajor axis (related to p+, 
the mean stellar densi ty), which can be d erived directly 
from the light curves (Soz zetti et al.ll2007f ). 

The stellar parameters were determined simultane- 
ously with the modeling of the light curves and radial 
velocities, as described next. We began by adopting the 
values of T c ff+, [Fe/H], and logg* from the SME analysis 
to fix the quadratic limb-darkening coefficients from the 
tabulations by Claret (2004), which are needed to model 
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the light curves (a,, &,). This modeling yields the proba- 
bility distribution of a/ R+ via a Monte Carlo approach, 
which is described fully in Section 13.31 We then used the 
a/ R* distribution together with Gaussian distributions 
for T c ff* and [Fe/H], with la uncertainties as reported 
previously, to estimate M* a nd R+ by compa rison with 
the stellar evolution models of lYi et all ([2001). This was 
done for each of 15,000 simulations. Parameter combina- 
tions corresponding to unphysical locations in the H-R 
diagram (4% of the trials) were ignored, and replaced 
with another randomly drawn parameter set. The re- 
sulting stellar parameters and their uncertainties were 
determined from the a posteriori distributions obtained 
in this way. 

In particular, the resulting surface gravity of log g* = 
4.34±0.03 is somewhat different from that derived in the 
initial SME analysis, which is not surprising in view of 
the possible strong correlations among T e fj*, [Fe/H], and 
the surface gravity. Therefore, in a second iteration we 
adopted this value of logg* and held it fixed in a new 
SME analysis, adjusting only T e ff*, [Fe/H], and wsini. 
This gave T eff * = 6158 ± 80 K, [Fe/H] = +0.17 ± 0.08, 
and wsini = 3.5 ± 0.5kms _1 , which we adopt as fi- 
nal atmospheric values for the star. Finally, we re- 
peated the calculation of the stellar mass and radius with 
these values and the stellar evolution models, and ob- 
tained M*=1.218± 0.039 M e , R*=l. 237 ± 0.054 R Q and 
i i =1.97±0.22LQ, along with other parameters summa- 
rized in Table [3] 

The model isochrones from lYi et al.l (|2001l) for metal- 
licity [Fe/H]=+0.17 are plotted in FigureHl with the final 
choice of effective temperature T e ^* and a/R+ marked, 
and encircled by the la and 2a confidence ellipsoids. 



TABLE 3 
Stellar parameters for HAT- P- 16 




6000 

Effective temperature [K] 



Fig. 4. — Stellar isochrones from lYi et al.l 11200 If ) for metallicity 
[Fc/H]=+0.17 and ages 0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0, 2.5 and 
3.0 Gyr (left to right). The final choice of T^ff* and a/R* are 
marked and encircled by the la and 2cr confidence ellipsoids. The 
initial value of T e ff + and a/R* from the first SME analysis and 
subsequent light curve analysis is marked with a triangle that is 
barely offset from the final value (indicated by a filled circle). 

The stellar evolution modeling provides color indices 
that may be compared against the measured values as 
a sanity check. The best available measurements are 
the ne ar-infrared magnitud es from the 2MASS Cata- 
logue (jSkrutskie et all l2006f ). J 2M ass = 9.850 ± 0.022, 
#2MASS = 9.623 ± 0.022 and K 2M ass = 9.553 ± 
0.020; which we have converted to the photometric sys- 
tem of the mod els (ESO) using the transformations by 
ICarpenterl (|2001l ). The resulting measured color index is 



Parameter 


Value 


Source 


T off * (K) 


6158 ± 80 


SME a 


[Fe/H] (dex) . . . 


+0.17 ±0.08 


SME 


vsini (kms ) 


3.5 ±0.5 


SME 


fmac (kms" 1 ) . 


4.61 


SME 


^mic (kms" 1 ).. 


0.85 


SME 


7RV (kms" 1 ) .. 


-16.83 ±0.19 


DS 




0.2166 


SME+Claret b 


bi 


0.3617 


SME+Claret 


M* (M ) 

R* («©) 

log 9* (cgs) .... 

i* (£©) 

V (mag) 


1.218 ±0.039 
1.237 ±0.054 
4.34 ±0.03 
1.97 ±0.22 
10.812 


YY+a/fl*+SME c 

YY+a/iJ*+SME 

YY+a/i^+SME 

YY+a/iJ*+SME 

TASS 


M v (mag) 


4.03 ±0.13 


YY+a/i?*+SME 


K (mag.ESO) 
M K (mag,ESO) 

Age (Gyr) 

Distance (pc) . . 


9.596 ±0.021 
2.74 ±0.10 
2.0 ±0.8 
235 ± 10 


2MASS+Carpenter d 
YY+a/iJ*+SME 
YY+a/i^+SME 
YY+a/_R*+SME 



a SME = "Spectroscopy Made Easy" package for analy- 
sis of high-resolution spectra Valcnti & Piskunov (1996). 
These parameters depend primarily on SME, with a small 
dependence on the iterative analysis incorporating the 
isochrone search and global modeling of the data, as de- 
scribed in the text. 

b SME+Claret = Based on the SME analysis and ta- 
blcs of quadratic limb-darkening coefficients from Claret 

JHDD- 

c YY+a/iJ*+SME = YY2 isochrones HYi et al.|[200lf) . 
a/Rir luminosity indicator, and SME results. 
d Based on the relations from Carpenter (2001). 

J-K = 0.319±0.032. This is within la of the predicted 
value from the isochrones of J—K = 0.33±0.02. The dis- 
tance to the object may be computed from the absolute 
K magnitude from the models (M K = 2.74 ± 0.10) and 
the 2MASS K s magnitude, which has the advantage of 
being less affected by extinction than optical magnitudes. 
The result is 235 ± 10 pc, where the uncertainty excludes 
possible systematics in the model isochrones that are dif- 
ficult to quantify. 

3.2. Excluding blend scenarios 

We performed a line bisector anal ysis of the observed 
spectra following iTorres et~a l. (2007) to explore the pos- 
sibility that the main cause of the radial velocity varia- 
tions is actually distortions in the spectral line profiles 
due to stellar activity or a nearby unresolved eclipsing bi- 
nary. The bisector analysis w as performed on the Keck 
spectra as described in § 5 of iBakos et a l. (2007a|); the 
spectra from the NOT were analyzed in a similar manner. 

The resulting bisector span variations are plotted in 
Figure [2] and have a low amplitude compared to the or- 
bital semi-amplitude and show no correlation with the 
radial velocities. We therefore conclude that the velocity 
variations are real, and that the star is transited by a 
close-in giant planet. 

3.3. Global modeling of the data 

Our model for t he follow-up light curve s used analytic 
formulae based on [Mandcl fc Agoll ([2002) for the eclipse 
of a star by a planet, where the stellar flux is described 
by quadratic limb-darkening. The limb darkening coef- 
ficients are based on the SME results ( Section [3Tj) and 
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the tables provided bv lClaretl (|2004| ) for the i band. The 
transit shape was parametrized by the normalized plan- 
etary radius R p /R+, the square of the impact parameter 
b 2 , and the reciprocal of the half mid-ingress to mid- 
egress duration of the transit Q/R*. We chose these 
parameters because of their simple geometric meanings 
and the fact that they show negligible correlations (see 
iBakoset al.ll2009D . Our model for th e HATNet data 
was th e simplified "P1P3" version of the Mandcl fc Agoll 
(2002) analytic functions (an expansion by Legendre 
polyn omials), for the reasons described in IBakos et al.l 
(2009). Following the formalism presented by Pal (2009), 
the RV curve was modeled as an eccentric Keplerian orbit 
with semiamplitude K, and Lagrangian orbital elements 
(k,h) — ex (cosw, sinw) 

We assumed that there is a strict periodicity in the 
individual transit times. In practice, we assigned the 
transit number N tr = to the first high quality follow-up 
light curve gathered on 2009 September 10. The adjusted 
parameters in the fit were the first transit center observed 
with HATNet, T Ci _284, and the last high quality transit 
center observed with the FLWO 1.2 m telescope, T c . + is. 
The transit center times for the intermediate transits 
were interpolated using these two epochs and the N t r 
transit number of the actual event (|Bakos et al.ll2007bh . 
The model for the RV data contains the ephemeris infor- 
mation through the T Ci _284 an d T c _ + i$ variables. Alto- 
gether, the 8 parameters describing the physical model 
were T Ci _ 2 84, ^c.+is, R P /R*, b 2 , (,/R*, K, k — ecosw, 
h = esinw. Nine additional parameters were related to 
the instrumental configuration. These are the three in- 
strumental blend factors i?i ns t,i of HATNet (one for each 
of the three fields), which account for possible dilution 
of the transit in the HATNet light curve due to the wide 
(20" wide FWHM) PSF and possible crowding, the three 
HATNet out-of-transit magnitudes, Mo,HN,i, an d three 
relative RV zero-points 7 ro ij (one each for the Keck, 
high- resolution FIES, and medium- resolution FIES ob- 
servations) . 

We extended our physical model with an instrumen- 
tal model that describes the systematic variations of the 
data. This was done in a simil ar fashion to the analysis 
presented in lBakos et ail (|2009| ). The HATNet photome- 
try has already been EPD- and TEA-corrected before the 
global modeling, so we only considered systematic correc- 
tion of the follow-up light curves. We chose the "ELTG" 
method, i.e. EPD was performed in "local" mode with 
EPD coefficients defined for each night, and TEA was 
performed in "global" mode using the same set of stars 
and TEA coefficients for all ni ghts. The under l ying p hys- 
ical model was based on the iMandel fc Agoll (|2002T ) an- 
alytic formulae, as described earlier. The five EPD pa- 
rameters were the hour angle (characterizing a monotonic 
trend that linearly changes over time) , the square of the 
hour angle, and the stellar profile parameters (equivalent 
to FWHM, elongation, position angle). The functional 
form of the above parameters contained six coefficients, 
including the auxiliary out-of-transit magnitude of the 
individual events. The EPD parameters were indepen- 
dent for all four nights, implying 24 additional coeffi- 
cients in the global fit. For the global TEA analysis we 
chose 20 template stars that had good quality measure- 
ments for all nights and on all frames, implying an addi- 



tional 20 parameters in the fit. Thus, the total number 
of fitted parameters was 17 (physical model) + 24 (local 
EPD) + 20 (global TFA) = 61. 



The joint fit was performed as described in lBakos et al.1 
(2009). We minimized \ 2 m the parameter space by us- 
ing a hybrid algorithm, combining the do wnhill simplex 
method (AMOEBA; see lPress et~aLlll992ft with the clas- 
sical linear least squares algorithm. Uncertainties for the 
parameters were derived usi ng the Markov Chain Monte- 
Carlo method (MCMC, see lFordll200fil IPa1l2009T ) . 

The eccentricity of the system appeared as significantly 
non-zero (k = -0.030 ± 0.003, h = -0.021 ± 0.006). 
The best-fit results for the relevant physical parame- 
ters are summarized in Table [4] We also list the RV 
"jitter," which is a component of assumed astrophysi- 
cal noise intrinsic to the star that we add in quadrature 
to the RV measurement uncertainties in order to have 
X 2 /dof = 1 from the RV data for the global fit. In addi- 
tion, some auxiliary parameters (not listed in the table) 
are: T c ,_ 2 84 = 2454297.5154 ± 0.0009 (BJD), T Cj+ i 8 = 
2455135.8554 ± 0.0003 (BJD), 7rcl = 4.5 ± 3.8ms -1 , 
7rd,FiEShi = -434.9 ± 4.2ms~ 1 (FIES, high resolution), 
7rci FiESmcd = —442.9 ± 2.8 m s _1 (FIES, medium resolu- 
tion), B mstrA23 = 0.77 ±0.08, B instrA63 = 0.93 ±0.04, 
Bi n str,i64 = 0.96 ± 0.04. Note that these gamma ve- 
locities do not correspond to the true center of mass 
RV of the system, but are only relative offsets. The 
true systemic velocity of the system is RV=— 16.83 ± 
0.19 km s^ 1 found by cross-correlating the observed spec- 
trum against a library template spectrum. The planetary 
parameters and their uncertainties can be derived by the 
direct combination of the a posteriori distributions of the 
light curve, RV and stellar parameters. We found that 
the planet is fairly massive with M p = 4.193 ± 0.094 Mj, 
and compact with radius of R p — 1.289±0.066 i?j, yield- 
ing a mean density of p p = 2.42 ± 0.35 g cm~ 3 . The final 
planetary parameters are summarized at the bottom of 
Table H 

4. DISCUSSION 

We present the discovery of HAT-P-16b with a period 
of P = 2.7760 d, a mass of M p = 4.19 ± 0.09 Mj and a 
radius of R p = 1.29±0.07 i?j. Currently, there are only a 
handful of planets residing near HAT-P- 16b in the m ass- 
radius diagra m. CoRoT-2b (lAlonso et al.l 120081 ) and 
CoRoT-6b (|Fridlund et al. 2010) have similar masses and 
periods, b ut both have ecce ntricities consist ent with zero. 
HD80 606 (|Naef et al.H200l and HD17156 dFischer et al.1 
120071 ) are also similar in mass, but both planets have 
long periods and very eccentric orbits. WASP-lOb (P = 
3.09d, e = 0.06, M p = 2.96 Mj,R p = 1.08 Rj, see 
iChristian et~aTl l2009t Uohnson et all [20091 ) is probably 
the transiting planet which most resembles HAT-P-16b 
with a similar period and eccentricity, albeit slightly 
lower mass and radius. At V = 10.8 mag, HAT-P-16 
is one of the brightest stars known to host a ~ 4 Mj 
planet. 

We ha ve compared HAT-P -16b to the theoretical mod- 
els from iFortney et al.l (|2007f ) by interpolating the mod- 
els to the solar-equivalent semi-major axis of a equ i V = 
0.0294 ± 0.0014 AU, the result of which can be seen over- 
plotted in Figure \5\ We find that the mass and radius 
of HAT-P-16b are quite consistent with the 1 Gyr model 
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TABLE 4 
Orbital and planetary parameters 



Parameter 



Value 



Light curve parameters 

P (days) 2.775960 ± 0.000003 

T c (BJD) 2455027.59293 ± 0.00031 

T w (days) a 0.1276 ± 0.0013 

T12 = T 34 (days) a 0.0150 ± 0.0014 

a/R* 7.17 ±0.28 

C/-R* 17.73 ±0.10 

Rp/R* 0.1071 ± 0.0014 

■ 2 nl qo+0.063 

uia ' : >-0.069 

b = acosi/R^ 0.439toog| 

i (deg) 86.6 ± 0^7 

RV parameters 

K (ms" 1 ) 531.1 ±2.8 

fc RV b -0.030 ± 0.003 

/i RV b -0.021 ± 0.006 

e 0.036 ± 0.004 

u 214 ±8° 

RV jitter (ms" 1 ) 8.0 

RV rms from fit (ms" 1 ) 10.0 

Secondary eclipse parameters 

T B (BJD) 2455028.929 ± 0.005 

T s 14 0.1234 ± 0.0020 

T s ,i2 0.0142 ± 0.0013 

Planetary parameters 

M p (Afj) 4.193 ± 0.094 

R p (iJj) 1.289 ± 0.066 

C(Mp,R p ) c 0.57 

p p (gem" 3 ) 2.42 ± 0.35 

a (AU) 0.0413 ± 0.0004 

log g p (cgs) 3.80 ± 0.04 

T cq (K) 1626 ± 40 

O d 0.220 ±0.011 

(F) (10 9 erg s" 1 cm" 2 ) c 1.58 ± 0.16 

a T14: total transit duration, time between first to last contact; 
T12 = T34: ingress/egress time, time between first and second, 
or third and fourth contact. 

Lagrangian orbital parameters derived from the global mod- 
eling, and primarily determined by the RV data. 
c Correlation coefficient between the planetary mass M p and 
radius R p . 

d The Safronov number is given by = h(Veac/V OI ^) 2 = 
(a/Rp)(A4 p /M ic ) fsee lHansen fc Barman]|2007l) . 
c Incoming flux per unit surface area. 

and conclude that HAT-P-16b is most likely a H/Hc- 
dominated planet. 

The radial velocity dataset used to characterize HAT- 
P-16b consist of three different velocity sets from two 
different telescopes (FIES medium and high resolution 
spectra from the NOT and HIRES spectra from Keck-I) , 
requiring additional free parameters in the global fit to 
account for the arbitrary velocity offset between the mea- 
surements. Nevertheless, the different observations yield 



very similar semi-amplitude and we find an rms from 
the best fit model of only 10.0 ms" 1 . The rich dataset 
consisting of 27 RV observations allows us to accurately 
determine the small non-zero eccentricity of the orbit as 
e = 0.036 with a high significance of 10c, with k 7^ at 
lOcr and h 5^ at 3cr. The eccentricity is consistent and 
clearly evident even when fitting separate orbits to the 
three different velocity datasets. 

Most short-period TEPs have eccentricities consis- 
tent with zero, as is expected because of circulariza- 
tion driven by tidal dissipation within the planet with 
a typical timescale sub stantiall y less than 1 Gyr (see 
I Jackson. Greenberg. fc Barned[2008l:lMardlindl2007[ ). In 
fact, the estimated circularization timescale for HAT-P- 
16b is r C i r ~ 0.03 Gyr (when the tidal dissipation pa- 
rameter is assumed to be Q p = 10 ), which is much less 
than the est imated ~ 2 Gyr a g e of t he system. On the 
other hand, iMatsumura et al.l (|2008T ) argue that circu- 
larization timescales could be significantly higher and 
the planets with eccentric or bits may simply be in the 
process of being circularized. lAdams fc Laughlml (2006) 
argue that for multiple systems with a hot Jupiter as 
the inner planet, secular excitation of the eccentricity by 
companion planets could explain the non-zero eccentric- 
ity of these systems. Thus the origin of the eccentricity 
of HAT-P-16b remains unclear. 

With a stellar rotation of 3.5 ± 0.5 kms^and an im- 
pact parameter of b ~ 0.43, HAT-P-16 is a favorable 
target for measuring the Rossiter-McLaughlin (RM) ef- 
fect. This effect can be used to determine the alignment 
of the planet's orbit al angular momen tum vector with 
the stellar spin axis (jWinn et al.l 120091 ) . The estimated 
amplitude of the RM effect is 47 ms -1 for HAT-P-16b. 
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